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Abstract 

The Monte Carlo method is the most popular technique to perform radiative transfer simulations in a general 3D geometry. The 
algorithms behind and acceleration techniques for Monte Carlo radiative transfer are discussed extensively in the literature, and 
many different Monte Carlo codes are publicly available. On the contrary, the design of a suite of components that can be used for 
the distribution of sources and sinks in radiative transfer codes has received very little attention. The availability of such models, 
with different degrees of complexity, has many benefits. For example, they can serve as toy models to test new physical ingredients, 
or as parameterised models for inverse radiative transfer fitting. For 3D Monte Carlo codes, this requires algorithms to efficiently 
generate random positions from 3D density distributions. 

We describe the design of a flexible suite of components for the Monte Carlo radiative transfer code SKIRT. The design is based 
on a combination of basic building blocks (which can be either analytical toy models or numerical models defined on grids or a 
set of particles) and the extensive use of decorators that combine and alter these building blocks to more complex structures. For a 
number of decorators, e.g. those that add spiral structure or dumpiness, we provide a detailed description of the algorithms that can 
be used to generate random positions. Advantages of this decorator-based design include code transparency, the avoidance of code 
duplication, and an increase in code maintainability. Moreover, since decorators can be chained without problems, very complex 
models can easily be constructed out of simple building blocks. Finally, based on a number of test simulations, we demonstrate 
that our design using customised random position generators is superior to a simpler design based on a generic black-box random 
position generator. 

Keywords: radiative transfer — methods: numerical — designing software — design patterns 


1. Introduction 

Radiative transfer, the process that describes radiation propa¬ 
gating through and interacting with matter, is a general problem 
that is encountered in all areas of astronomy (and far beyond). 
Due to its high dimensionality and nonlocal and nonlinear be¬ 
haviour, it is generally considered as one of the most challeng¬ 
ing problems in numerical astrophysics. 

Recent years have seen an impressive advancement of three- 
dimensional (3D) radiative transfer studies, thanks to the in¬ 
crease in computational power, the availability of a wealth of 
new observational constraints and the development of new al¬ 
gorithms. Among the different approaches available to solve the 
radiative transfer problem, the Monte Carlo method is generally 
the most popular one. The first Monte Carlo radiative transfer 
codes were developed more than four decades ago (e.g., Mattila 
1970; Roark et al. 1974; Witt 1977), and since then the method 
has steadily increased its market share in all fields where 3D ra¬ 
diative transfer is important. Many powerful Monte Carlo codes 
are available in various fields of computational astrophysics, in¬ 
cluding dust radiative transfer (Steinacker et al. 2013, and ref¬ 
erences therein), Lya line transfer (Verhamme et al. 2006; Di- 
jkstra et al. 2006; Laursen et al. 2009), ionising radiation trans¬ 
port (Ciardi et al. 2001; Maselli et al. 2003; Wood et al. 2004), 


neutron transport (Romano & Forget 2013) and neutrino radia¬ 
tion transport (Abdikamalov et al. 2012). These and many other 
papers and monographs discuss at length the main algorithms 
behind Monte Carlo radiative transfer and the various proposed 
improvements and acceleration techniques. 

A very different aspect of Monte Carlo radiative transfer 
codes that has received very little attention in the literature, is 
the setup of a suite of components that describe the distribution 
of the sources and sinks in the radiative transfer code (i.e., the 
objects that add radiation to or remove radiation from the ra¬ 
diation field). This is an aspect outside the real core of Monte 
Carlo radiative transfer problem. Virtually all of the available 
codes require the user to hard-code the model makeup for each 
distinct problem (although the results of hydrodynamical sim¬ 
ulations can generally be processed out of the box, given the 
appropriate import module). We argue, however, that it is very 
useful for a radiative transfer code to provide a suite of input 
models, or geometries, as built-in components. Such toy mod¬ 
els can provide a low-threshold introduction to new users. They 
also have a crucial role in a more scientific way: toy input mod¬ 
els are used to benchmark different codes (Ivezic et al. 1997; 
Pascucci et al. 2004; Pinte et al. 2009), to investigate the ef¬ 
fects of dust attenuation on the apparent structural parameters 
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of galaxies (Mollenhoff et al. 2006; Gadotti et al. 2010; Pas- 
trav et al. 2013a,b), or to investigate the physical impact of e.g. 
dumpiness or spiral arms (Witt & Gordon 1996; Bianchi et al. 
2000; Misiriotis et al. 2000; Semionov et al. 2006). Finally, 
they are essential for so-called inverse radiative transfer, i.e. 
when a parameterised radiative transfer model is fitted to obser¬ 
vational data (Xilouris et al. 1999; Bianchi 2007; Schechtman- 
Rook et al. 2012; De Geyter et al. 2013, 2014). 

In order to be useful for these goals, the suite of input mod¬ 
els should be diverse enough and contain models with different 
degrees of complexity, ranging from the metaphorical spherical 
cows to more realistic toy models that consider, for example, 
spheroidal and triaxial geometries, including bars, spiral arms 
and clumpy distributions. Setting up such a suite is more com¬ 
plex than it might appear at first sight. While each model is in 
principle completely determined by the 3D density distribution 
p(jc), the setup requires more than just an implementation for 
this density function. 

A crucial aspect of Monte Carlo radiative transfer codes is 
the emission of a multitude of simulated photon packets from 
random locations sampled from the source density distribution. 
So each model in the suite should contain a routine that gener¬ 
ates random numbers according to p(jc). Moreover, this random 
position generation needs to be very efficient, since the random 
positions are sampled extremely often. In this aspect, Monte 
Carlo radiative transfer simulations differ from other codes where 
random positions need to be sampled from arbitrary density dis¬ 
tributions, e.g. to generate the initial conditions for N-body or 
hydrodynamical simulations. The efficient generation of ran¬ 
dom positions from complex 3D density distribution is not a 
trivial task. 

Another aspect that needs consideration is the organisation 
of such a suite of input models. One could provide a param¬ 
eterised model that can cover every possible option, but this 
would quickly lead to an explosion of options that is hard to 
overview and maintain, and would inevitably contain substan¬ 
tial code duplication. 

In this paper, we describe how these issues are dealt with 
in the publicly available 3D Monte Carlo dust radiative transfer 
code SKIRT^ (Baes et al. 2003, 2011; Camps & Baes 2015). 
While we summarise the relevant information for the discus¬ 
sion here, an in-depth description of the architecture and over¬ 
all design of the SKIRT code can be found in Camps & Baes 
(2015). It should also be noted that the presented design issues, 
while discussed in the context of SKIRT and thus dust radia¬ 
tive transfer, are fully applicable to other Monte Carlo transport 
problems. 

This paper is laid out as follows. In Section 2 we review the 
general techniques to generate random vectors from multivari¬ 
ate distributions that are available in the specialised numerical 
analysis literature, but less so in the astrophysics community. In 
Section 3 we present the general setup of the suite of models for 
the SKIRT code. In Section 4 we describe the various building 
blocks currently present in the SKIRT code, and in Section 5 we 


present a number of decorators that can be used to combine and 
add complexity to simple building blocks. Some decorators are 
analysed in more detail, focusing on the methods used to effi¬ 
ciently generate random positions. In Section 6 we discuss the 
advantages of the decorator-style design of our suite of com¬ 
ponents, and we critically investigate an alternative option in 
which random positions are generated using a generic routine 
rather than customised generators. We show that, while such 
an approach is simpler to implement and might hence seem an 
attractive alternative, it cannot compete with our approach, in 
terms of accuracy and efficiency. Section 7 sums up. 


2. Multivariate random number generation techniques 

The generation of random numbers from univariate distribu¬ 
tions is a well-known topic in numerical analysis. A number 
of standard methods are widely used and clearly described in 
standard textbooks (e.g.. Press et al. 2007). Additional tech¬ 
niques that are less known in the (astro)physics community in¬ 
clude the acceptance-complement method (Kronmal & Peter¬ 
son 1981, 1984) and the Forsythe-von Neumann method (von 
Neumann 1951; Forsythe 1972). For an extensive overview 
of random number generation from univariate distribution, see 
Devroye (1986). Unfortunately, the generation of random vec¬ 
tors from multivariate distributions is much more complex. The 
only multivariate distributions from which random vectors can 
be generated directly are those where the density can be written 
as a product of independent univariate density distributions. In 
general, however, more advanced techniques are necessary. 


2.7. The inversion method 


The inversion method, also known as inversion sampling, is the 
most popular method for univariate generation problems. The 
basis of this method is the following: if y is distributed accord¬ 
ing to a density^ g(y), then the variable x defined as the solution 
of the equation y = y(x) is distributed according to the density 


fix) = g(y(Ji:)) 


ax 


( 1 ) 


If we now want to generate a number from a given density /(x), 
we can take a uniform distribution 


giy) = 



if 0 < y < 1, 
else. 


( 2 ) 


and set y(x) = F(x) where F is the cumulative distribution 
corresponding to the density /(x). The inversion method can 
hence be used to generate random variates with an arbitrary 
density, provided that the inverse function 7^“^(y) to the cumu¬ 
lative distribution F(x) is explicitly known. Classical examples 
include the exponential distribution, the Cauchy distribution, 
the Rayleigh distribution and the logistic distribution. 


^ SKIRT documentation: http://www.skirt.ugent.be 
SKIRT code repository: https://github.com/skirt/skirt 


^Throughout this paper, all densities, univariate or multivariate, are assumed 
to be properly normalised. 
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Formula (1) can directly be expanded to multiple dimen¬ 
sions: if y is distributed according to the joint probability dis¬ 
tribution g(y), then the vector jc defined as the solution of the 
vector equation y = j(jc) is distributed according to the joint 
distribution 


fix) = 



( 3 ) 


where \dy/dx\ is the Jacobian determinant. Probably the most 
famous application of this formula is the Box-Muller method to 
generate random normally distributed deviates (Box & Muller 
1958; Bell 1968). Another interesting special application is the 
case of a linear transformation. In this case, the transformation 
y = y(x) can be written as a matrix multiplication y = lix, and 
the distribution of x is 


f(x) = \li\g(lix) (4) 

with I H I the absolute value of the determinant of H. This 
kind of transformation is particularly useful for the generation 
of random vectors with a given dependence structure, as mea¬ 
sured by the covariance matrix (Scheuer & Stoller 1962; Barr 
& Slezak 1972). 

In general, however, it is not straightforward to use this 
identity (3) to construct a method that can be used to generate 
random positions from an arbitrary multidimensional distribu¬ 
tion. 


increases rapidly when going from one to more dimensions, 
which decreases the efficiency of the method. Moreover, the 
design of a suitable reference distribution becomes much more 
complicated. 

2.3. The composition method 

The composition method or probability mixing method is an¬ 
other important technique in both univariate and multivariate 
random number generation (Marsaglia & Bray 1964; Hermann 
et al. 2004). Rather than a method on itself to generate non- 
uniform random numbers, it is a principle to facilitate and speed 
up other random number generating methods. The simple idea 
behind composition is to decompose the distribution f{x) as a 
weighted sum, 

K 

fix) = Tj ^ifiix) ( 5 ) 

i=l 

with all fi normalised densities, and the weights Wi a proba¬ 
bility vector (i.e., all wi ^ 0 and 2/^/ = !)• To generate a 
random x from the distribution /(jc), we first generate a ran¬ 
dom integer number k from the discrete probability vector w/, 
and subsequently generate a random x from the density fkix). 
A prerequisite for this method to work is that, obviously, the 
decomposition can be done efficiently, and that the complexity 
of the problem is reduced by the decomposition. 


2.2. The rejection method 

Apart from the inversion method, the rejection method, also 
known as the acceptance-rejection method, is the most popular 
method to sample nonuniform random numbers from univariate 
distributions (von Neumann 1951; Devroye 1986). The basic 
idea behind the method is that, if one wants to sample a ran¬ 
dom number from a function /(v), one can sample uniformly 
from the 2D region under the graph f{x). More concretely, as¬ 
sume that f{x) ^ c/ref(x), where /ref(x) is another distribution 
from which random numbers are easily generated, and c is the 
so-called rejection constant (it obviously satisfies the condition 
c ^ 1). One then generates a uniform deviate ^ and a random 
number x from the reference distribution /ref(x), and calculates 
the quantity t = ^c/i-ef(x)//(x). This procedure is repeated until 
^ ^ 1, in which case x is the desired random number. 

One of the advantages of the rejection method is that it does 
not require that the cumulative distribution function be analyt¬ 
ically known, let alone be invertible. However, its effective¬ 
ness depends on how accurate / is approximated from above 
by c/ref. Less accurate approximation leads to a greater chance 
of rejection; on average c iterations of the loop are required 
before one successful random number is generated. Moreover, 
the reference distribution should be such that random numbers 
can be easily generated from it and that the computation of 
is simple. For a range of standard distributions, 
such as the gamma distribution and the Poisson distribution, ef¬ 
ficient reference functions can be constructed. 

The rejection method is by no means limited to univari¬ 
ate distributions, and can immediately be applied to multivari¬ 
ate generation problems. However, the rejection rate typically 


2.4. The conditional distribution method 

Finally, a powerful technique that applies only to multivariate 
random number generation is the so-called conditional distribu¬ 
tion method (Devroye 1986; Hormann et al. 2004). It is based 
on the Bayesian identity 

fiX\,X2)= fi{Xi)f2iX2\Xi) (6) 

which expresses that a joint distribution of two variables can be 
calculated as the product of the marginal distribution of the for¬ 
mer variable and the conditional distribution of the latter. Ex¬ 
panding it to multiple dimensions, 

/(X) = /l(Xi)/2(X2|Xi)---/Ar(XA^|Xi,...,XA^_i) (7) 

The beauty of this technique is that it reduces the multidimen¬ 
sional generation problem to a sequence of independent uni¬ 
variate generation problems. The main drawback is that it can 
only be used when much detailed information is known about 
the distribution /(x). In particular, it is by far not always the 
case that marginal and conditional distributions are easily cal¬ 
culated in closed form. The standard textbook example of this 
technique is the multivariate Cauchy distribution. 

3. General setup of the SKIRT Geometry suite 

SKIRT is a multi-purpose 3D Monte Carlo dust radiative trans¬ 
fer code that is mainly used to simulate dusty galaxies and ac¬ 
tive galactic nuclei (e.g., Baes & Dejonghe 2002; Baes et al. 
2010; Stalevski et al. 2012; De Looze et al. 2012, 2014). The 
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code was designed as a highly modular software, with a par¬ 
ticular consideration for the development of a flexible and easy 
user interface and the use of proven software design principles 
as described by Camps & Baes (2015). 

The SKIRT code offers a wealth of configurable features 
that are ready to use without any programming. In particular, 
the SKIRT code is equipped with a suite of input model com¬ 
ponents, the so-called Geometry classes, that can be used to 
represent distributions of either sources or sinks. Essentially, 
the suite consists of a number of classes that all inherit from 
an abstract Geometry class, for which the C-r-r class interface 
looks like 

class Geometry { 
public: 

Geometry 0; 
virtual ~Geometry(); 

virtual double density(Position x) const = 0; 
virtual Position generatePositionO const = 0; 

} 

Both the density and generatePosition functions are pure 
virtual functions, which means that each model in the suite 
should provide a routine that implements the (normalised) den¬ 
sity p{x), and a routine that generates random positions accord¬ 
ing to this density. For example, the interface of a concrete 
GonereteGeometry class that contains only a single parame¬ 
ter p would look as follows 

class ConcreteGeometry : public Geometry { 
public: 

ConcreteGeometry(double p) {_p=p} 
double density(Position x) const; 

Position generatePositionO const; 
private: 
double _p; 

} 

For the design of this suite of input models, a naive option 
would be to provide a set of parameterised models that con¬ 
tain free parameters to control every possible option. A typical 
standard component of such a suite would be a very generalised 
Plummer model, with free parameters setting the location of 
the centre, the orientation with respect to the coordinate sys¬ 
tem, the scales, the flattening describing potential triaxiality, 
possible degrees of dumpiness, etc. This approach has a num¬ 
ber of strong disadvantages. Clearly, it would quickly lead to 
an explosion of different options that is hard to overview. It 
would inevitably lead to substantial code duplication (nearly all 
the code would need to copied if we would consider a Hern- 
quist profile instead of a Plummer profile) and code that is vir¬ 
tually impossible to maintain (the code for all models would 
need to be updated if a new feature is added or altered). Also 
concerning efficiency, such a design would not be optimal. In¬ 
deed, a random position generating routine for such a gener¬ 
alised Plummer model is hard to construct and is certainly much 
less efficient than the simple routine that is possible for a plain 
spherical Plummer model (using the inversion method). 


To overcome these problems, we adopted a completely dif¬ 
ferent approach that is much simpler but still provides the flex¬ 
ibility and functionality needed to set up complex models. This 
is achieved using a combination of simple base models on the 
one hand, and so-called decorator geometries on the other hand. 
Decorator geometries are not real models on their own, but they 
apply modifications upon other models in interesting ways, fol¬ 
lowing the Decorator design pattern. In object-oriented pro¬ 
gramming, a decorator attaches additional responsibilities to an 
object dynamically and provides a flexible and powerful alter¬ 
native to subclassing for extending functionality (Gamma et al. 
1994). 

In our present context, a decorator geometry is a special 
kind of geometrical model (i.e., it is also a C-r-r class in the gen¬ 
eral Geometry suite) that takes one or more other components 
and adds a layer of complexity upon them. Simple examples 
of decorators that can easily be implemented are the relocation 
of the centre of a given component, or the rotation of a given 
model with respect to the coordinate system. More complex 
decorators deform a spherical model to a triaxial one, or add a 
spiral perturbation to an axially symmetric model. The advan¬ 
tage of the decorator approach is clear: each decorator needs 
to be implemented only once, and can then be applied to any 
possible model. 

In the following two sections we describe the building blocks 
in the SKIRT Geometry suite, and a number of decorator ge¬ 
ometries that can alter these building blocks to more complex 
structures. 

4. The SKIRT Geometry building blocks 

The SKIRT Geometry suite contains a limited number of ele¬ 
mentary input models, which can be used either as elementary 
toy models, or as building blocks for more complex geometries. 
For each of them, the density can be expressed as a simple ana¬ 
lytical function, and the generation of random positions reduces 
to three independent univariate generation problems. 

Apart from these analytical components, SKIRT also offers 
the possibility to set up components in which the geometry of 
sources and/or sinks is defined by means of particles or on a 
grid. In particular, SKIRT can import a snapshot from a (mag- 
neto)hydrodynamical simulation. One obvious goal would be 
to post-process a hydrodynamical simulation in order to cal¬ 
culate the observable multi-wavelength properties of the simu¬ 
lated objects (e.g., Jonsson et al. 2010; Schartmann et al. 2014). 
In this case, it would be sufficient to just read in the geometry of 
the snapshot as it is, and start the radiative transfer simulation. 
More generally, however, it would be useful if these numerical 
components would be at a similar level as the analytical com¬ 
ponents. This would open up the possibility to decorate them 
and combine them with other analytical and/or numerical com¬ 
ponents to more complex models. 

4.1. Analytical components 

A first group consists of spherically symmetric models. In spher¬ 
ical symmetry, the generation of a random position from the 
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3D density p{x) = p(r) simplifies to the generation of a ran¬ 
dom azimuth from a uniform distribution, a random polar angle 
from a sine distribution, and a random radius from the uni¬ 
variate density f{r) = 47ir^p(r). The SKIRT suite includes 
the most popular models used to represent shells, star clus¬ 
ters, early-type galaxies, or galaxy bulges, such as a power-law 
shell model (Ivezic et al. 1997), the Plummer model (Plum¬ 
mer 1911), the y-model (Dehnen 1993; Tremaine et al. 1994), 
the Sersic model (Sersic 1963; Ciotti & Berlin 1999; Baes & 
Gentile 2011), and the Einasto model (Einasto 1965; Retana- 
Montenegro et al. 2012). Within this family of models, many 
other famous profiles are contained, including the Hernquist 
model (Hernquist 1990), the Jaffe model (Jaffe 1983) and the 
7^1/4-model (de Vaucouleurs 1948). 

A second group of elementary input models consists of ax- 
isymmetric models in which the density is a separable function. 
The standard example with a density distribution separable in 
cylindrical coordinates is the double-exponential disc that is the 
de facto standard to represent disk galaxies in radiative transfer 
studies (Xilouris et al. 1997; Baes et al. 2003; Bianchi 2007). 
How random positions can be drawn from this distribution is 
discussed in Appendix A of Baes et al. (2003). An example of 
an axisymmetric model where the density is separable in spher¬ 
ical coordinates is the torus model that is often used to represent 
the dusty tori around stars or active galactic nuclei (Collison & 
Eix 1991; Granato & Danese 1994; Stalevski et al. 2012). 

4.2. Components based on smoothed particles 

The first group of numerical components in SKIRT are defined 
as a set of smoothed particles. This is mainly useful when 
we want to use the output of a smoothed particle hydrodynam¬ 
ics (Lucy 1977; Monaghan 1992; Springel 2010b) simulation. 
In spite of claims that the technique suffers from fundamen¬ 
tal problems (Agertz et al. 2007; Bauer & Springel 2012), it 
is still the most popular hydrodynamics technique, especially 
for cosmological simulations of galaxy formation (e.g., Guedes 
et al. 2011; Eeldmann & Mayer 2015; Schaye et al. 2015). The 
output of an SPH snapshot consists of a set of “particles” (or 
rather anchor points in a co-moving grid), each characterised 
by a large suite of physical quantities. 

As far as SKIRT is concerned, most of these physical quan¬ 
tities are irrelevant. An SPHGeometry component in SKIRT 
is essentially defined by a list of N smoothed particles and the 
assumed smoothing kernel W{r, h), with each smoothed parti¬ 
cle characterised by a position Xj, a fractional mass mj and a 
smoothing length hj. The total density at an arbitrary position 
X is then given by 

N 

p(x) = J^mjW(\x-xj\,hj) (8) 

j=i 

In practice, the kernels used in SPH simulations almost always 
have a finite support (e.g., Monaghan & Lattanzio 1985; Des- 
brun & Gascuel 1996; Muller et al. 2003), and in this case only 
a relatively small number of terms in the sum have a non-zero 
contribution. 


The SPHGeometry class in SKIRT employs smoothing ker¬ 
nel implementations optimised for each specific task. The ge¬ 
ometry’s density at a given location is calculated according to 
equation (8) using a finite-support cubic spline kernel, so that 
the operation can be limited to particles that potentially over¬ 
lap the location of interest. To facilitate this process, the setup 
phase of the simulation places a rough grid over the spatial do¬ 
main and constructs a list of overlapping particles for each grid 
cell. As described in section 3.3 of Camps & Baes (2015), a 
further optimisation is provided to calculate the mass within a 
given box (a cuboid lined up with the coordinate axes), as an al¬ 
ternative to sampling the density in various locations across the 
volume of the box. In this case, the calculation uses the analyti¬ 
cal properties of a scaled and clipped Gaussian kernel, designed 
to approximate the cubic spline kernel, to directly determine the 
mass in the box. This optimisation accelerates the density cal¬ 
culation for typical cartesian grids by an order of magnitude. 

On the other hand, the generation of random positions sam¬ 
pled from the geometry’s density distribution is rather straight¬ 
forward, thanks to the composition method. The first step is 
the choice of a random smoothed particle, based on a discrete 
distribution where each particle is weighted by its relative mass 
contribution. The second step is generating a random position 
according to the distribution of the chosen particle. The cur¬ 
rent implementation samples a random number from a Gaus¬ 
sian smoothing kernel with infinite support using the inversion 
method. 

In the future we may provide a suite of smoothing kernels, 
allowing the user to select the kernel that best fits the SPH snap¬ 
shot being imported. The methods used for calculating the den¬ 
sity and generating a random radius would obviously depend 
on the actual shape of the kernel. 

In addition to what is described above, SKIRT has facili¬ 
ties to associate a spectral energy distribution with each SPH 
star particle based on properties such as age and metallicity, or 
to associate a dust mass with each SPH gas particle based on 
properties such as temperature and metallicity. In the former 
case, the geometry will weigh the particles by luminosity, for 
each wavelength separately, and in the latter case, the geometry 
will weigh the particles by dust mass. A detailed description of 
these features is beyond the scope of this paper. 

4.3. Components based on hierarchical grids 

Apart from smoothed particle hydrodynamics, the main other 
technique that is used to perform hydrodynamics simulations is 
Eulerian mesh-based hydrodynamics (Stone & Norman 1992; 
Eryxell et al. 2000). A fundamental ingredient of this technique 
is the use of grids with adaptive mesh refinement. Eulerian 
AMR simulations are used in virtually all fields of astrophysics, 
and Monte Carlo codes have been adapted to work directly on 
these hierarchical grids (Kurosawa & Hillier 2001; Niccolini & 
Alcolea 2006; Saftly et al. 2013, 2014). 

The HierarchicalGridGeometry in SKIRT reads in snap¬ 
shots defined by density fields discretised on hierarchical carte¬ 
sian grids and converts them to the format of the other com¬ 
ponents in the SKIRT Geometry suite. Calculating the den- 
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sity at an arbitrary position conies down to identifying the cell 
that contains this position and returning the density associated 
with that cell. As hierarchical grids have the structure of a tree, 
isolating the correct cell is a straightforward and computation¬ 
ally cheap operation. Generating random positions from such a 
component is similar to the case of the smoothed particles, and 
is based on the composition method. We first generate a random 
cell from a discrete distribution where each cell is weighted by 
its relative mass. Secondly, we determine a random position 
within the chosen cell; as the cells in a hierarchical cartesian 
grids are cuboids, this is a trivial task. 

4.4. Components based on Voronoi grids 

Recently, a new Lagrangian technique that solves the hydro¬ 
dynamics equations on a moving, unstructured Voronoi grid is 
gaining popularity. It is claimed to avoid some of the difficulties 
of smoothed particle hydrodynamics on the one hand and Eu- 
lerian grid-based hydrodynamics on the other hand. This tech¬ 
nique has been used for many years in the computational fiuid 
dynamics community (Mavriplis 1997), and a number of novel 
codes based on this principle have recently been developed in 
the astrophysics community (Springel 2010a; Duff ell & Mac- 
Fadyen 2011). Moving mesh hydrodynamics is mainly applied 
to simulations of galaxy formation and evolution (e.g., Mari- 
nacci et al. 2014; Vogelsberger et al. 2014). 

SKIRT contains a VoronoiGridGeometry class that con¬ 
verts a snapshot from a Voronoi hydrodynamical simulation (or 
any other density field defined on a Voronoi grid) to a regular 
SKIRT Geometry component. Due to the nature of a Voronoi 
grid, the only necessary input is the list of all the generating 
sites and the associated densities; it is hence not necessary to 
store all the vertices and sides of each of the cells. Based on the 
generating sites, SKIRT constructs the corresponding unique 
Voronoi grid using the public Voro-r-r library (Ry croft 2009). 

The density routine essentially works in the same way as 
for the case of the hierarchical grids: it comes down to identi¬ 
fying the cell that contains the given position and returning the 
density associated to this cell. In the case of a Voronoi grid, 
however, the cell identification is not as simple as in a hierar¬ 
chical tree structure of cartesian grid cells. Due to the nature 
of Voronoi grids, this is essentially a nearest neighbour search. 
Rather than looping over all possible sites, SKIRT implements 
an approach using cuboidal blocks, as explained in detail in 
Camps et al. (2013). This task could be optimised even further 
using more advanced techniques based on space partitioning 
structures such as k-d trees or R-trees (Friedman et al. 1977; 
Guttman 1984; Flaw et al. 2010). 

Also the generation of random positions works essentially 
in the same way as for hierarchical grids. The first step is iden¬ 
tical: we generate a random cell from a discrete distribution 
where each cell is weighted by its relative mass contribution. 
The second step, generating a random position from within the 
chosen cell, is significantly more complex than in the case of a 
cuboidal cell. To the best of our knowledge, there are no ded¬ 
icated techniques to generate a random point from a Voronoi 
cell. There are two possible options. 


Table 1: An overview of the geometry decorators currently implemented in 
SKIRT, referencing the section in which they are presented. Geometry decora¬ 
tors can be applied to basic geometry building blocks, and can be chained and 
combined to create complex structures. 


5.1 

Offset 

applies an arbitrary offset to any geometry 

5.2 

Rotate 

applies an arbitrary rotation to any geometry 

5.3 

SpheCavity 

Crop 

carves out a spherical cavity from any geometry 
crops any geometry to a given box; a variation on 
the Cavity decorator 

5.4 

Combine 

combines two geometries into a single geometry 

5.5 

Triaxial 

Spheroidal 

transforms any spherical geometry to a triaxial form 
transforms any spherical geometry to a spheroidal 
form; a special case of the Triaxial decorator 

5.6 

Spiral 

applies spiral arm structure to any axisymmetric ge¬ 
ometry 

5.7 

Clumpy 

replaces a portion of the mass in any geometry by 
randomly placed clumps 

6.2 

Foam 

replaces the model-specific random position genera¬ 
tor by a generic routine based on the Foam library 


The first option is to partition the cell into a set of tetrahe- 
dra, subsequently select a random tetrahedron from a discrete 
distribution where every tetrahedron is weighted by its relative 
volume, and finally generate a random position from the se¬ 
lected tetrahedron. Specific algorithms are available for both 
the tetrahedrisation of convex polyhedra (Edelsbrunner et al. 
1990; Max 2002) and the generation of random positions from 
a tetrahedron (Rocchini & Cignoni 2000). 

The second option, which is more simple and which we 
have adopted in SKIRT, is to use the rejection technique. As 
the reference distribution we use a uniform density in a cuboidal 
volume, defined as the 3D bounding box of the cell. As Voronoi 
cells are convex polyhedra, this bounding box is directly ob¬ 
tained when the vortices of the cell are known (these are calcu¬ 
lated using the Voronoi grid setup). Extensive testing has shown 
that, depending on the distribution of the generating Voronoi 
sites, the average ratio of the volume of the bounding box of a 
Voronoi cell over the actual cell volume is about 3 to 4. This 
ratio immediately represents the average rejection rate for the 
random position generation. 

5. The SKIRT Geometry decorators 

In this section we describe a number of decorator geometries 
that can be applied on the building blocks described in the pre¬ 
vious section in order to convert them to more complex struc¬ 
tures; see Table 1 for an overview. The implementation of the 
density of a decorator geometry is usually not a major problem; 
the main challenge is to implement the routine that generates 
random positions from a decorator geometry, so this is what we 
focus on. 

5.1. Offsets 

The Of f setGeometryDecorator decorator in SKIRT applies 
an arbitrary offset a to any density distribution. If the original 
density is Ps(j), the new density is simply p(x) = p^(x - a). 
Generating random positions is equally simple: we generate a 
random y from the original density distribution and return x = 
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y + a. The C++ implementation of the class in SKIRT looks 
like: 

class OffsetGeometryDecorator : public Geometry { 
public: 

OffsetGeometryDecorator(Geometry* g, Position a) 
{_g = g; _a = a;} 
double density(Position x) const 
{return _g->density(x-a);} 

Position generatePositionO const 
{return _g->generatePosition()+a;} 
private: 

Geometry* _g; 

Position _a; 

} 

5.2. Rotation 

Similarly, the RotateGeometryDecorator decorator rotates 
any density distribution. If the rotation is characterised by the 
orthonormal matrix H, the new density is p{x) = ps(Hjc). To 
generate a random position from this new density, we generate a 
random y from the original density and rotate it over the inverse 
rotation matrix, i.e. x = H^j. 

5.3. Cavities 


5.5. Triaxial geometries 

As a first more complex decorator, we consider the case of a 
triaxial decorator geometry, which converts a spherically sym¬ 
metric density distribution into one that is stratified on concen¬ 
tric and confocal ellipsoids. More concretely, assume that we 
have a spherically symmetric density distribution Ps(j) = PsCv), 
we then consider its triaxial counterpart 


Pix) = — Ps 
pq 



( 11 ) 


Such triaxial models are discussed and used extensively to de¬ 
scribe the stellar distribution of elliptical galaxies and galactic 
nuclei and the shape of galactic haloes (e.g.. Stark 1977; Mer¬ 
ritt & Fridman 1996; Trujillo et al. 2002; van den Bosch et al. 
2008). Oblate and prolate spheroidal distributions are just a 
special case of these triaxial models in which p = \. 

It is clear that the density (11) cannot be written as a product 
of independent univariate density distributions. We can use the 
conditional distribution method. In order to develop a general 
recipe for generating random positions, we start by rewriting 
the probability distribution according to (11) in spherical coor¬ 
dinates. 


A third simple decorator, the CavityGeometryDecorator, car¬ 
ves out a cavity from another density Ps(j). In formula form, we 
have 

i O jc in cavity 

PsW (9) 

with;^ the fraction of the mass of the original density located in 
the cavity. This decorator is useful to represent density distribu¬ 
tions of dust close to a star or active galactic nucleus, where the 
dust has been cleared due to sublimation. To generate random 
positions from this new density distribution, we just generate a 
random position from the original density distribution and re¬ 
ject is when it is located in the cavity. This is, in fact, an almost 
trivial application of the rejection technique, where the origi¬ 
nal density assumes the role of the reference function, and the 
rejection constant isc = 1/(1-;^). 


5.4. Composition 

The CombineGeometryDecorator combines two density dis¬ 
tributions into a single distribution according to 

Wipi{x) + W2p2{x) 

p{x) = - (10) 

Wi + 1+2 

withpi,p 2 the original distributions and i+i, 1+2 their respective 
weights in the composite distribution. Generating random posi¬ 
tions for this new density distribution is a trivial application of 
the composition method. 


r^sinO 

f{r,e,(p) = -ps 

pq 


sin^0 (sin^^ + cos^^) cos^0 


( 12 ) 

According to formula (7), we now rewrite this expression as 

/(r,^,0) = /(0)/(^|0)/(r|^,0) (13) 


After some calculation, one finds for the marginal distribution 
for 0 the simple expression 


m = 


j_ p 

'2n sitV<p + p2 cos^<p 


(14) 


Note that this marginal distribution is independent of the spe¬ 
cific shape of the density profile and only depends on the batten¬ 
ing parameter p (and forp= 1, it reduces to a simple uniform 
distribution, as expected for a spheroidal distribution). Gener¬ 
ating a random azimuth from this density can be done using the 
standard inversion method; the corresponding cumulative dis¬ 
tribution is 


F((f)) = arctan 

ZTT 


tan0 

P 


(15) 


and is readily inverted. Once this random 0 has been deter¬ 
mined, we can generate a random polar angle 6 from the condi¬ 
tional distribution 


f(0\(p) 


sin 6 I sin^O 
2s \ s^ 


+ cos 6 


-3/2 


(16) 
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where we have set s^ = q^(sin^(p + p^ cos^0)//7^. Also this dis¬ 
tribution is independent of the specific shape of the density pro¬ 
file and only depends on the battening parameters p and q (and 
















the already randomly determined azimuth). Again, generating 
a random polar angle from this distribution can be done using 
the standard inversion method: the corresponding cumulative 
distribution is 




Vl + 


(17) 


and is again readily inverted. Finally, we consider the condi¬ 
tional distribution for r. 


fir I 0) = 


/ sin^^ 


.3/2 


1 Ps 

V 


\ 

If we apply a linear transformation from r to 


-r cos'^6 


r sin 6 


y= -■ 
q 


+ cos^O 


(18) 


(19) 


where (R,(p,z) are the usual cylindrical coordinates, and the 
function ^(R, 0) is defined as 


fiR, (p) = Cn sin^^ 


n ( IniR/Rp) 
2 \ tan p 



(23b) 


The factor is n normalisation constant that guarantees that 
the density is normalised, and is given by 




V^r(A-r 1) 

m+k) 


(24) 


This general spiral arm perturbation is a generalisation of the 
models used by Misiriotis et al. (2000) and Schechtman-Rook 
et al. (2012), and allows for an arbitrary number of arms n, pitch 
angle p, spiral perturbation weight w, and arm-interarm size 
ratio (controlled by the parameter N). For A = 1 the density 
reduces to a more simple perturbation with equal arm-interarm 
size ratio. 


the corresponding density for m is simply 

f(y \0,tp) = 47ry VsCv) (20) 

which is nothing but the distribution for a random radius from 
the original density Ps(j). The last step in the generation of a 
random position is to generate a random position y from the 
original spherical density distribution Ps(j), and transform the 
radius y = Ijl of this position to a new radius r using the linear 
transformation (19), in which we use the previously determined 
values for 6 and 0. 

While this method based on conditional probabilities has a 
certain beauty, there is, actually, a simpler and more efficient 
method that is based on the inversion method. Indeed, consider 
the simple transformation 

(y\,y2,y3) = (xux2lp,x3iq) (2i) 


or in matrix format 


y = lix 


(1 0 

0 2 

p 

[o 0 


O' 

0 

1 


( 22 ) 


p(x) =psiR,z){l + wsin 


(ln(R/Ro) 

' - 

\ tan p 


-4> 


(25) 


For the general 3D density distribution (23), we can again use 
the conditional distribution method, and write 


/(R,0,z) = /(R,z)/(0|R,z) (26) 

Since the spiral perturbation is such that it averages out along 
every single circle around the z-axis, i.e., 

^(R,0)d0=l (27) 



we trivially find 


f(R,z) = 27iRp,(R,z) (28) 

This is exactly the 2D probability density distribution that cor¬ 
responds to the density of the axisymmetric density distribution. 
We hence simply generate a random position x from our orig¬ 
inal density distribution and save the radius R and the height 
z. In order to generate a random azimuth 0, we consider the 
conditional distribution /(01 R, z), for which we find directly 


According to the law of linear transformations (4), if y is dis¬ 
tributed according to the spherical density p^iy), x will be dis¬ 
tributed exactly according to the triaxial density (11). An easy 
way to generate random positions from a triaxial counterpart 
of a spherical symmetric density is thus simply generating a 
random position y from the spherical density distribution, from 
which the desired position can easily be calculated as jc = H^j. 

5.6. Spiral arm structure 

As a second complex decorator, we consider a logarithmic spi¬ 
ral arm perturbation that can be applied to any axisymmetric 
density distribution, 

p{x) = Ps(/?, z) [(1 - w) + w^{R, (p)] (23a) 


mR,z) 


(1 - w) -I- wfiR, 0) 
In 


(29) 


This univariate density is independent of the shape of the origi¬ 
nal axisymmetric density profile; it only depends on the param¬ 
eters of the spiral perturbation and randomly determined radius 
R. For the general expression (23b), the standard inversion tech¬ 
nique is not easily applicable; the cumulative distribution corre¬ 
sponding to the density (29) can be calculated analytically, but 
the resulting expression cannot be inverted analytically. One 
could resort to a numerical inversion, but a better approach is 
to use the rejection technique. We can use a simple uniform 
density as the reference distribution, /ref(0) = 1/2;:. The rejec¬ 
tion constant c should ideally be chosen as the smallest value 
for which /(01R, z) ^ <^/ref(0) for all values of 0. Since the 
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maximum value of the perturbation function 0) is Ca^, we 
find 

c=\ + {Cn-\)w (30) 

Actually, as an alternative to the conditional probability tech¬ 
nique, we could also simply have used the rejection technique 
from the beginning. The 3D density p{x) satisfies the condition 
p(jc) ^ cps(jc), so if we simply use the original density Ps(jr) as 
our reference distribution, we have essentially the same algo¬ 
rithm with the same rejection constant. The difference between 
this 3D rejection technique and the previous version, where we 
used the rejection technique only for the conditional distribu¬ 
tion for the azimuth, is a matter of efficiency. In the former 
version only one random position needs to be generated from 
the original density, and it the rejection is only to determine 
the azimuth. In the latter version, an entirely new position is 
generated for every rejection, which is less efficient. 

5.7. Clumpy models 

As a last example, we consider the ClumpyGeometryDecorator 
class that turns any density distribution Ps(j) into a clumpy ana¬ 
logue, i.e., a density distribution in which a fraction of the mass 
is distributed “smoothly” as in the original density distribution, 
and the remaining fraction is locked up into compact clumps. 
The density of this distribution can be written as 

f ^ 

p(x) = (1 - /ci)Ps(x) + h) (31) 

where /d is the fraction of the mass in clumps, N is the number 
of clumps, and the positions Xt are the central positions of the N 
clumps, determined randomly according to the original density 
function ps. The function W(r, h) is the smoothing function that 
sets the distribution of matter within every single clump, with 
h a characteristic length scale. This can in principle be any 
spherical 3D distribution; in practice we use either a uniform 
density sphere or one of the compact support kernels that are 
used in smoothed particle hydrodynamics simulations. 

It is straightforward to generate random positions from this 
distribution. It is a direct application of the composition method, 
and the case is very similar to smoothed particle building blocks 
discussed in Section 4.2. 

6. Discussion 

6.1. Advantages of the decorator-style design 

As indicated in Section 3, a major advantage of the decorator- 
style approach is that each decorator needs to be implemented 
only once, and can then be applied to different models. This 
avoids the need to implement complex geometries with many 
different possible features and heavy code duplication. That 
alone is already an argument strong enough to justify its use. 

Another strong advantage is the flexibility offered by this 
approach. The decorators discussed in Section 5 transform an 
input model into a more complex model. Such a decorated 
model is often still relatively simple: it typically adds one layer 


of complexity onto the model that is being decorated (think of 
a rotated exponential disc or a triaxial Sersic model). If we 
want to generate toy models that can be used as idealised rep¬ 
resentations of, for example, a spiral galaxy, various layers of 
complexity have to be added. The decorator-style approach is 
ideal for this purpose: decorators can be applied also to models 
that have already been decorated. In other words, decorators 
can easily be chained or nested. 

The left-hand panel of Figure 1 shows a projected image of 
a relatively complex toy model, constructed using the Geometry 
suite in SKIRT. The model consists of a bulge component and 
a clumpy spiral disc. Each of these two components have a 
very simple spherical y-model as their starting point, on which 
a chain of decorators was applied. The former was first deco¬ 
rated into a triaxial model, the latter was turned into an axisym- 
metric disc by applying a very strong flattening, then a spiral 
perturbation was added, and a fraction of its mass was turned 
into clumps. 

In decoration chains the successive decorators do not nec¬ 
essarily need to be different. It is also possible to chain the 
same decorator several times, i.e. to apply them recursively. 

A powerful example of this nesting is a repeated application 
of the ClumpyGeometryDecorator. The right-hand panel of 
Figure 1 shows an application of this principle on a simple tri¬ 
axial Plummer model. In this case, we have applied four nested 
applications of the decorator, in which the number of clumps in¬ 
creases in each level, whereas the smoothing length decreases. 
This algorithm results in a fractal density distribution that is 
self-similar over an order of magnitude in scale (Elmegreen 
1997; Mathis et al. 2002; Indebetouw et al. 2006). 

The possibility to chain decorators, including a repeated ap¬ 
plication of the same decorator, facilitates the construction of 
very complex models out of simple building blocks. Analyti¬ 
cal components can easily be combined with numerical compo¬ 
nents based on smoothed particles or grids, and decorators can 
be applied regardless of the underlying component type. It is, 
for example, possible to add a smooth triaxial bulge to an irreg¬ 
ular disc structure defined using SPH particles, or to carve out 
a cavity from a complex hydrodynamics system and locate a 
small nuclear structure there to simulate the effects of an AGN 
on the large-scale structure of a galaxy. 

6.2. The use of a generic random position generator 

One could argue that the setup that we have chosen is unnec¬ 
essarily complex. An alternative approach would be a design 
where we only provide the density for each component in the 
suite, and where the generation of random positions is executed 
by a generic routine rather than a geometry-specific routine for 
each component. Such a suite could still use the decorator de¬ 
sign pattern, but would be simpler to implement. The main 
challenge is the design of a generic routine that can gener¬ 
ate random positions corresponding to an arbitrary 3D density 
function. 

This is in principle possible: there are a few multivariate 
black-box random number generation libraries available, in¬ 
cluding Foam (Jadach 2000,2003) and RanLip (Beliakov 2005a,b), 
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Figure 1: Two examples of the chained application of decorators. The left panel is a simulated image of a two-component model consisting of a triaxial bulge, 
and a flattened disc model on which a three-armed logarithmic spiral perturbation and a clumping decorator were applied. For both the bulge and the disc 
component, the starting model was a simple spherical 7 model. The model to the right is a triaxial Plummer model, with four recursive applications of the 
ClumpyGeometryDecorator decorator. The result is a model with a fractal structure. In both models 10^ photon packages were released to create the images, and 
the noise in the images is negligible (i.e., all the features in the image are real). 



run time (s) 

Figure 2: Results of the timing tests for the test models #1 through #5, using the 
model-speciflc and the generic Foam-based random position generators with 
25,000 through 400,000 cells. The horizontal bars show the simulation run¬ 
time in seconds, split into the fraction devoted to the setup of the grid (red) and 
the actual time spent on random position generation (green). 


based on the so-called composition-rejection method, a com¬ 
bination of the composition and rejection methods. In a first 
phase, the exploration phase, the domain of a distribution is 
partitioned into different cells. In the generation phase, the re¬ 
jection technique is used on each cell, where typically a con¬ 
stant density function is used as the reference distribution. The 
advantages of this approach are that, in principle, it can be ap¬ 
plied for all distributions and in any dimension, and that only 
the density of each component is required. The main disadvan¬ 
tage is the complexity and the overhead, both in memory and 
run-time, that are linked to the exploration phase. 

In order to test whether or not our design choice of provid¬ 
ing a customised random position generator for each concrete 
Geometry class (and for each decorator, in particular) is justi¬ 
fied, we have set up a comparison between our routines and a 
generic black-box routine. We have implemented a new decora¬ 
tor class, FoamGeometryDecorator, that replaces the model- 
specific random position generator by a generic routine based 
on the Foam library. Foam is a self-adapting cellular code that 
uses iterative binary splitting, using either simplexes or hyper¬ 
rectangles, to subdivide the configuration space. 

We set up a simple suite of test models, for which we gen¬ 
erate random positions using both the model-specific generator 
and the generic Foam generator. Our suite consists of a se¬ 
quence of models with an increasing level of complexity. It es¬ 
sentially follows the different steps in building the model shown 
in the left panel of Figure 1 . The following models are consid¬ 
ered 

1. a spherical 7 = ^ model, 

2. model #1 flattened to an axisymmetric disc, 

3. model #2 with spiral perturbation, 

4. model #3 with a fraction turned into clumps, 

5. model #4 with a triaxial 7 = ^ bulge added. 
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For each model, we generate 10^ random positions (we do not 
use any photon propagation or detection as would be used in 
Monte Carlo radiative transfer simulation in order to isolate the 
random position generation process). The timings were done 
using a single core and with averages over multiple test runs. 
One crucial aspect for the Foam-based generator is the choice 
of the Foam parameters, in particular of the number of cells 
in the grid. A larger number of cells implies a higher compu¬ 
tational cost of the exploration phase on the one hand, but a 
better approximation of the density and smaller rejection rates. 
The ideal choice of this parameter is impossible to determine in 
a general way. We consider 5 different values for the number 
of cells in the Foam for every model, ranging from 25,000 to 
400,000. 

The results are illustrated in Fig. 2, where, for each model, 
we give the total run time and the contributions for the setup 
of the foam and the actual time spent on generating the random 
positions. 

A first major conclusion that can be drawn from these re¬ 
sults is the efficiency of the customised random position gen¬ 
eration routines. From a simple spherical model, to a complex 
model composed of two components that both have been deco¬ 
rated multiple times, the computational cost increases only by 
a factor two. Moving from a spherical to a triaxial model does 
not affect the efficiency of the generation of random positions at 
all (not a surprise, as it just implies two simple multiplications, 
as shown in Section 5.5). The biggest jump in timing occurs 
when spiral structure is added to the model (i.e., from model #2 
to #3), because the random position generation is based on the 
rejection method, which has more overhead compared to the 
inversion and composition methods. The addition of clumping 
(from model #3 to #4) hardly implies an increase in computa¬ 
tion time. Note that the addition of a bulge even decreases the 
computation time. In model #5 half of the positions are gen¬ 
erated according to the density of the relatively simple triaxial 
bulge, whereas in model #4 they all are generated according to 
the more complex clumpy spiral structure. 

A second major conclusion is that the generic Foam-based 
random position generator is significantly less efficient than the 
customised generator. For the most simple models (#1 and #2), 
the difference in speed is a factor 3-8, depending on the number 
of cells used. The increase in total run time for the Foam-based 
models with increasing number of cells is mainly due to the 
time necessary to set up the foam. The time necessary to gen¬ 
erate the random positions also increases, but not as strongly. 

For the more complex models, the efficiency of the Foam- 
based generator decreases even more. The biggest jump in tim¬ 
ing now occurs when we add clumping to the model, because 
the calculation of the density for the clumpy models is com¬ 
putationally very demanding (each density evaluation requires 
a sum over a large number of clumps). For models #4 and 5, 
the Foam-based generators are at least an order of magnitude 
slower than the customised routines. Interestingly, the run time 
of these simulations does not increase with the number of cells. 
For simulations with few cells, the setup phase is relatively 
quick, but the generation phase is very inefficient. The reason 
is that the cells are relatively large with a strong variation in 


density values, and hence the rejection rates are uncomfortably 
large. If the number of cells in the grid increases, the setup time 
naturally increases, but the actual position generation becomes 
more efficient because the average rejection rates are smaller. 

Moreover, the inefficiency of the Foam-based generator is 
not the only problem, but also accuracy is an issue, in particu¬ 
lar for complex models with multiple local maxima and strong 
gradients, such as models #4 and #5. If the number of cells is 
limited and the cells hence rather large, both the total mass and 
the local maximum within a cell are very hard to guess (they are 
determined by random sampling the cell). As a result, both the 
relative weights of the different cells and the rejection constants 
within the individual cells are poorly determined. This will lead 
to an artificial smoothing and wrong results. Figure 3 shows 
a sequence of images corresponding to model #4, created us¬ 
ing the customised random position generator (top left) and us¬ 
ing the Foam-based one with different values for the number of 
cells (remaining panels). The customised generator generates 
an image that reveals all the details that are expected, includ¬ 
ing the sharp density gradients and the contrast between arm 
and interarm regions. The only limitation in the image is the fi¬ 
nite resolution due to the pixel size and the unavoidable Poisson 
noise. The Foam images on the other hand show a clear signa¬ 
ture of degradation that gradually decreases when the number 
of cells grows. For the grids with 25,000 and 50,000 cells, and 
even for the one with 100,000 cells, the effects of the grid are 
clearly visible, and the individual clumps and the spiral struc¬ 
ture are insufficiently resolved. For the grid with 400,000 cells, 
the individual clumps are well resolved, but the sharpness of 
the spiral arms is still under-resolved, especially in the central 
regions. 

Since both the accuracy and the efficiency of the generic 
Foam-based random position generator cannot compete with 
the customised versions, we conclude that it is worth investing 
in the latter, and that our design choice is justified. 

7. Conclusions 

We have described the design of a suite of components that can 
be used to model the distribution of sources and sinks in the 
Monte Carlo radiative transfer code SKIRT. Our main conclu¬ 
sions are the following: 

• The availability of a well-designed suite of input models, 
with enough variety and different degrees of complex¬ 
ity, in a publicly available Monte Carlo code has a strong 
added value. Such models can serve as toy models to 
test new physical ingredients, or as parameterised mod¬ 
els for inverse radiative transfer fitting. They also provide 
a low-threshold introduction to new users, since models 
with differing degrees of complexity can be run without 
any coding at all. 

• Each model is in principle completely determined by the 
3D density distribution p{x). In order to be used in Monte 
Carlo radiative transfer, however, each model should con¬ 
tain a routine that generates random numbers according 
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Figure 3: A comparison of the Monte Carlo generated images for test model #4, based on 10^ photon packages. For the top-left panel, we have used the customised 
random position generator based on the chain of decorators. For the other panels, we have used the generic Foam-based generator, with the number of cells in each 
panel indicated in the top-left corner. 


to this p{x). Finding the most suitable algorithm to im¬ 
plement this is the most challenging aspect of the design 
of a suite of input models. 

• The design of the SKIRT Geometry suite is based on 
a combination of basic building blocks and the exten¬ 
sive use of decorators. The building blocks can either 
be simple analytical components, or they can be numer¬ 
ical components defined as a set of smoothed particles 
or on a hierarchical or unstructured Voronoi grid. The 
Geometry decorators combine and alter these building 
blocks to more complex structures. 

• The different multivariate random number generation tech¬ 
niques that exist in the specialised numerical analysis lit¬ 
erature can be used to efficiently implement complex dec¬ 
orators, for example those that add triaxiality, spiral struc¬ 
ture or dumpiness to other models. Decorators can be 
chained without problems and essentially without limita¬ 
tion. Different layers of complexity can hence be added 
one by one. The result is that very complex models can 
easily be constructed out of simple building blocks, with¬ 
out any coding at all. 


• From the software design point of view, decorators have 
many advantages, including code transparency, the avoid¬ 
ance of code duplication, and an increase in code main¬ 
tainability. This is a clear example that adhering to proven 
software design principles pays off, even for small and 
mid-sized projects. 

• Finally, we demonstrate that our design using customised 
random position generators is superior to a simpler suite 
design based on a generic black-box random position gen¬ 
erator. Using a suite of test simulations with increasing 
complexity we demonstrate that our customised random 
number generators are more accurate and more efficient. 
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